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We present preliminary development of an approach for simulating high Reynolds num- 
ber steady compressible flow in two space dimensions using a Cartesian cut-cell finite 
volume method. We consider both laminar and turbulent flow with both low and high 
cell Reynolds numbers near the wall. The approach solves the full Navier-Stokes equations 
in all cells, and uses a wall model to address the resolution requirements near boundaries 
and to mitigate mesh irregularities in cut cells. We present a quadratic wall model for 
low cell Reynolds numbers. At high cell Reynolds numbers, the quadratic is replaced with 
a newly developed analytic wall model stemming from solution of a limiting form of the 
Spalart-Allmaras turbulence model which features a forward evaluation for flow velocity 
and exactly matches characteristics of the SA turbulence model in the field. We develop 
multigrid operators which attain convergence rates similar to inviscid multigrid. Inves- 
tigations focus on preliminary verification and validation of the method. Flows over flat 
plates and compressible airfoils show good agreement with both theoretical results and 
experimental data. Mesh convergence studies on sub- and transonic airfoil flows show con- 
vergence of surface pressures with wall spacings as large as ~0.1% chord. With the current 
analytic wall model, one or two additional refinements near the wall are required to obtain 
mesh converged values of skin friction. 


I. Introduction 

C ARTESIAN embedded-boundary grids have proven extremely useful for inviscid flow simulations around 
complex geometries. Their primary strengths include their accuracy, rapid turnaround time and level of 
automation. The use of cut cells at the intersection of the grid and the geometry has been well-studied, and 
the numerical issues of discretizations, stiffness and convergence for these irregular cells are well understood. 

These issues are less well understood for high Reynolds number flows on cut-cell meshes. Since the 
Navier-Stokes equations involve one higher derivative, the numerical issues are more delicate. The literature 
shows difficulties extracting smoothly varying quantities such as skin friction due to irregularity of the cut 
cells, along with loss of accuracy and numerical stiffness. In addition, the resolution requirements needed to 
compute aerodynamic flows at Reynolds numbers of interest are daunting. The inability of Cartesian meshes 
to (easily) refine anisotropically in the wall normal direction is clearly a challenge for these methods. 

There are several approaches found in the literature for dealing with the mesh irregularity and resolution 
requirements of Cartesian meshes. The most obvious is to avoid cut cells altogether, and use layers of 
conformal cells before switching to a background Cartesian mesh. 1-3 Some authors accomplish this through 
generation of a body-fitted grid in the near-body region and then changing to a background Cartesian grid. 
Alternatively one can start with a cut-cell mesh, remove the layers of cells adjacent to the body, and then 
drop normals to the body to create the body-aligned mesh layers . 4,5 Both methods rely on body-fitted 
near- wall layers to avoid the large cell counts associated with isotropically refined Cartesian grids, but give 
up the simplicity and robustness of a pure Cartesian approach. An alternative idea is to use the Cartesian 
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mesh down to the wall. This can be done using either a finite volume 6,7 or finite difference 8,9 scheme, or 
other integration methods 10,11 but all will have to contend with the non-aligned boundary. Some methods 
use an immersed boundary 12, 13 instead of explicitly cutting the cells that intersect geometry. For any of 
these approaches to be ultimately practical, some additional technique must be used to alleviate the cell 
count problem of isotropic refinement for high Reynolds number flows. 

Each of these approaches has benefits and drawbacks. We have chosen to use Cartesian cut cells so 
that all mesh faces remain Cartesian, placing emphasis on the ability to automatically generate meshes 
around complex configurations. In fully Cartesian meshes, the volume mesh is decoupled from the surface 
triangulation. This simplicity is one of the great strengths of embedded-boundary methods and is key for 
automation. Using only Cartesian faces also allows highly optimized flow solvers. 

Recent papers 14-20 present promising results using wall function or PDE-based wall model approaches 
to mitigate the inefficiency of isotropically refined Cartesian meshes in a boundary layer. This is also the 
approach we are taking in our current work, suitably modified for finite volume meshes with Cartesian 
cells explicitly cut by the embedded boundary. In this preliminary work, we model near-wall behavior with 
either a simple quadratic representation of the solution (low cell Reynolds numbers) or with an analytic wall 
function (high cell Reynolds numbers) as our initial wall model for turbulent flow. While wall functions are 
sufficient for this initial exploratory work, ultimately we envision developing a more complete PDE-based 
subgrid wall model. 

In this paper we lay the groundwork in several steps. We first develop a Navier-Stokes solver for laminar 
flows, by developing an accurate discretization for the cut cells and mesh interfaces. This is discussed in 
Section II and Appendix V. It entails an understanding of grid irregularity at the cut cells for second- 
derivative terms. Of particular importance is the use of a quadratic to compute the values in the cut 
cells (for both the solution and its derivatives) needed for the flux computations in laminar flow. This 
work builds upon our earlier efforts 21,22 where the governing equations are integrated to steady-state using 
a second-order cell-centered finite volume scheme. Section II also presents the RANS equations and the 
Spalart-Allmaras turbulence model. The coupling of the wall model to the finite volume scheme used in the 
turbulent flow computations is also presented here. Section III discusses fortification of the basic multigrid 
solver for viscous flow. Section IV presents two-dimensional computational examples of both laminar and 
turbulent flow, including flat-plate boundary layers on both coordinate-aligned and non-aligned meshes, and 
compressible airfoils. 


II. Discretization of the Navier-Stokes Equations 


In integral form in two space dimensions the compressible Navier-Stokes equations can be written 
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where the vector Vg is proportional to the gradient of the temperature \7q = — /cVT, the Prandtl number 

is given by Pr = fic p /k , and p, is computed using Sutherland’s law. 23 We take Pr = .72. 

We first present the basic finite volume scheme used on the regular Cartesian cells that make up most 
of the volume mesh. Only the viscous terms will be discussed, since the inviscid terms have been previously 
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described. 22 Following that we present the modifications needed for the viscous discretization in the irregular 
cells. Of particular importance here will be the method used to compute the gradient at the faces of these 
cut cells. 

A. Cartesian Cell Discretization G Q 

\ U y, \ U V;, , i 

The discretization of the viscous terms on un- 
cut Cartesian cells is straightforward. At the 
midpoint of each Cartesian face, the gradient of 
both velocity components is needed. Consider 
for simplicity an x face. We ensure stability 
and accuracy by taking u x « Ut+1, ^~ Uz,J at that 
face. This avoids the pitfall of simply averaging 
the left and right cell-centered gradients, (al- 
ready computed using second-order central dif- 
ferences for the inviscid terms in a second-order 
finite volume scheme). This would lead to an 
unstable 5-point approximation to the second 
derivative with a decoupled stencil. The trans- 
verse derivative however, in this case u y at the x face, can stably be taken as the average from the adjacent 
cells, u y « .5 {u Vi j + u Vi+1 d ). This discretization is commonly found in the literature for regular Cartesian 
meshes. 24,25 

At mesh refinement interfaces, a least-squares gradient computation that uses only face neighbors is 
linearly exact, and so approximates the solution to second order. We use this gradient to recenter the 
solution to compute the fluxes at mesh interfaces in the same way described in (21) for the cut cells. Since it 
is not centered, this difference approximation gives only a first-order accurate approximation to the derivative 
at an interface. The transverse component of the gradient is still taken to be an average of the transverse 
gradients on either side of the face, this time weighted by distance from the face. 

B. Cut-Cell Discretization for Navier-Stokes 

At cut faces, a recentering procedure described 
in Figure 2 is used to reconstruct the velocity 
from the cell centroid to the perpendicular line 
through the face centroid. Appendix A presents 
an accuracy study of this discretization com- 
pared with two other popular methods using a 
Poisson equation model problem. This study 
shows second-order accuracy in both the L\ and 
max norms, which is on a par with the best of 
these methods. 

When the face lies between two cut cells, 
both cells follow this procedure and the deriva- 
tive can be computed. If one of the cells is un- 
cut, the face itself must not be cut, and only one 
cut cell needs to recenter the solution, since the 
full Cartesian cell centroid and the face centroid 
must already be aligned. Pressure at the wall is 
obtained by simple reconstruction from the cell centroid to the surface centroid of the wall in each cut cell. 

Mesh irregularity through the cut cells introduces non-smooth truncation errors along walls. Accordingly, 
special procedures must be devised for accurate reconstructions of velocity derivatives at wall boundaries. For 
meshes with low cell Reynolds numbers at the wall, we obtain accuracy by using a higher-order (quadratic) 
polynomial basis for reconstruction of the data in the wall- normal direction. In RANS simulations (discussed 
later) with high cell Reynolds numbers, this quadratic basis is replaced by an appropriate wall model. 



Figure 2. Recentering of the solution from cell centroids 
(A and C) to points (B and D) on the line perpendicular 
to the midpoint of the cut-cell edge. A simple difference 
(B-D)/Ax approximates the derivative. 
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Figure 1. Illustration of stable viscous flux computation at 
face between cells i and i + 1. 
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1. Mesh Irregularity at the Cut Cells 

Mesh irregularity adjacent to geometric 
boundaries has been a major challenge for 
an accurate discretization of the viscous 
terms in the Navier- Stokes equations us- 
ing Cartesian mesh methods. Since the 
viscous terms include derivative quanti- 
ties, irregularity affects both the com- 
puted viscous fluxes and the skin friction 
computed at the wall. Accurate skin fric- 
tion, in particular, has been a challenge 
and has stymied many viscous Cartesian 
efforts. Figure 3 displays the problem 
and is representative of many of the re- 
sults found in the literature. 26 28 This 
figure shows the skin friction Cf along 
a ReL = 5000 flat plate rotated 15° to 
the Cartesian mesh at a free stream Mach 
number of 0.5. The magenta line indicates the Blasius solution, and the symbols represent the data extracted 
from cut cells along the plate’s surface. Skin friction data using the viscous discretization method with lin- 
ear reconstruction outlined in the preceding section have been colored by number of neighbors. Like similar 
examples in the literature, the skin friction data appear quite noisy around a mean which roughly follows 
the Blasius result for Cf. In this particular example, however, the sorting of data by number of neighbors 
reveals a stratification of the noise with cell type. On a 15° flat plate, 2-neighbor cells are always the smallest 
cut cells while 4-neighbor cells are always the largest. As a result, the patterns of the symbols associated 
with 2-, 3- and 4- neighbor cells graphically illustrate the link between mesh and stencil irregularity and skin 
friction noise. 

Figure 4 reveals the source of this irregularity. The figure shows a non-aligned body with cut cells a, 
5, c and d. To the right, we sketch a nonlinear profile U(y) marked with the cell-centroid data for cells a 
through d and slopes U' (a) through U'(d). With a non-linear profile, the slopes at a-d are obviously not 
(■ dU/dy) wa u , and the “noise” staircasing through the profile in Fig. 3 is not surprising. In the skin friction 
plot, the cell-centroid gradients are being taken along the wall, even though the centroids in this non-uniform 
region are at different distances from the wall. 

This observation suggests an obvious path forward: given cell centroid data at the first and second 
cells off the wall, reconstruct a quadratic function through these data, and evaluate the wall slope using 
this reconstruction. Indeed, this reconstruction can provide all relevant slopes needed by the discretization 
stencil outlined in the preceding section - the wall slopes as well as gradients at cell centroids. 



Figure 3. Noise in the skin friction profile along surface of the 
plate when using piecewise linear gradients in the cut cells, as in 
the literature. 26-28 



Figure 4. Illustration of a quadratic interpolant through the cut cells on a non-aligned mesh. 
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2. Near- Wall Quadratic Reconstruction for Low Cell Reynolds Numbers 

As suggested above, it is a simple matter to fit a quadratic to three values, particularly since both the u and 
v velocity is zero at the wall. Using the notation of Figure 4, we construct a line through cell l’s centroid, 
normal to the patch of the wall contained in cell 1. Cell 2 is the neighbor across the face intersected by 
this normal. Its centroid is a distance y 2 in the normal direction from the point 0 at the wall. Strictly 
speaking, data in cell 2 should be recentered tangentially as well. However, since streamwise gradients in 
high Reynolds number flows are generally several orders of magnitude smaller than wall-normal gradients, 
the current implementation uses data at 2 directly. 

The resulting formula is 

U(y) = yA + y(y-y 1 )B (3) 

with 

u 2 -ui _ 

A=^, and B = (4) 

y 1 V2 yo 

U (y) can be differentiated in the wall normal direction giving 

U'(y) = A + (2y- yi )B. (5) 

This gives a slope at the wall of 

(fO ~ V(0) = A — yiB (6) 

\ d Vj wall 

We recenter the u and v velocities independently. The slope at the wall provides the wall shear stress 
for both viscous fluxes and output skin friction, after rotating into the tangential direction. In addition, 
the gradient of the quadratic evaluated at the cell centroid replaces the least-squares gradients previously 
computed and is used in both inviscid and viscous flux balances. The quadratic reconstruction is essentially 
compensating for the stencil irregularity near the wall, and it is used both in the cut cells and also in the 
first layer of interior cells (cell 2 in Fig. 4). Although these cells are full hexahedra, they also have boundary- 
dependent stencil irregularity since they are adjacent to the cut cells. The use of quadratic reconstruction 
mitigates the mesh irregularity while maintaining strict conservation. 2 * * 5 

Of possible concern in using the quadratic to discretize wall shear is that it may adversely impact the 
allowable time step. Appendix B contains a stability analysis for a model ID heat equation to determine if a 
quadratic boundary condition reduces the CFL limit. Results using GKS stability analysis 29 show that using 
forward Euler in time and central differencing in space results in a CFL limit of 0.43. The standard scheme 
using linear reconstructions is stable with a CFL limit of 0.5, so the impact on the timestep is minimal. Our 
choice of time step at the cut cells is sufficiently conservative that this reduction is not an issue. 

C. RANS equations 

Compressible high Reynolds number turbulent flows are modeled using the Favre- averaged Navier-Stokes 
equations (herein still referred to as the RANS equations). They are of the same form as (1) except for the 
stresses, which are empirically modeled in turbulent flow. Here we use the Boussinesq assumption relating 
the Reynolds stress tensor to known flow properties. The net effect of this is the addition of an eddy- viscosity 
parameter / i t and a turbulent Prandtl number Pr t which we take to be 0.9, giving shear stresses 

2 

T XX ^ (/^ T Vy) 

T X y = (l* + lh)(u y + V x ) = T yx ( 7 ) 

2 

T yy = g (^ + Lt)(^ v y ~ u x) 

and turbulent heat flux Vg = (^ + ^^)VT. 

We model fit using the Spalart-Allmaras turbulence model 30 in fully turbulent form (i.e. no transition 
or tripping terms. We use what is referred to as the baseline model described in Oliver’s 2008 dissertation. 31 

2 One area for potential improvement is that the quadratic is fit to the cell-averaged value as if it were the pointwise solution, 

as is common for second-order finite volume schemes. Strictly speaking a higher order scheme should account for this, but this 

error is small. 
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Added to this are terms for the so-called negative model , which modifies some of the functional forms in the 
production and destruction terms for negative values of v to help push the turbulent viscosity back to the 
positive region. 31,32 

The resulting equation in the incompressible formulation is 
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where we made the usual assumption that ^ is small and can be neglected in rearranging the diffusion 
terms as above. For this somewhat modified negative model we take 
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and S = ? where Xl is the magnitude of vorticity. For the destruction term we have 
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with fi t = pi)f v i and the usual definitions of f w , g, r and y, and the other constants. 

This negative model improves on and expands the robustness measures in the original and baseline 
models. 30,31 Oliver applied the negative model in a GMRES framework. 31 In our multigrid context, we 
found some of the modifications for negative v overly steep and reduced them somewhat while still following 
the same design principles. Reducing them substantially improved the depth of convergence of our multigrid 
smoother. We have not done extensive study in arriving at the final functional forms. All equations - Navier 
Stokes plus turbulence model - are solved in a fully coupled manner as described in Section III. 

One additional change was necessary that is not commonly discussed in the literature. The advective 
terms in the turbulence model are usually implemented using a first-order scheme to help preserve positivity. 
This was too diffusive on non-streamwise-aligned meshes. The first-order scheme was replaced by a second- 
order method with linear reconstruction to discretize the advective terms. 

The Spalart-Allmaras turbulence model uses the normal distance to the wall in the destruction term (the 
d in 11)) for each cell in the mesh. We use Sethian’s fast marching method 33 to compute this additional 
quantity. This is mostly standard, with only a few modifications needed for our multi-level cut-cell mesh. 
This is incorporated into our mesh generator as a preprocessing step before computing the flow. 


1. Cut-Cell Wall Model for High Cell Reynolds Numbers 

One of the obvious challenges facing Cartesian methods for RANS simulations are the untenable cell counts 
associated with resolving the viscous stresses all the way down to the wall. Subgrid wall modeling is based 
on the observation that very near the wall, momentum transfer is governed by viscous stresses and the 
near wall flow is distinct from the outer flow in which inertial forces dominate. Wall functions, which are 
essentially algebraic wall models, are based on this observation. While wall- function development has been 
active 34 since the 80’s, recent work has focused on extending them for numerical approaches using subgrid 
modeling for both body-fitted 15-17 and non- body- fitted grids. 18,19 An advantage of these approaches is 
that the streamwise flow gradients can be taken from the underlying grid cells, thus removing many of the 
restrictions of algebraic wall models with built-in assumptions about the underlying flow conditions (zero 
pressure gradient etc.) These two-layer models can also couple more flexibly to the outer RANS flow. This is 
done for example in the diffusion model of Bond and Blottner, 20,35 where a system of odes is solved instead 
of assuming an analytic form for the wall function. The issues then become how far from the wall does the 
coupling occur, and what equations are solved in the wall layer. 
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While our eventual goal is to develop a complete subgrid model, in this paper we first examine the 
issues that arise when extending the use of wall functions to non-coordinate-aligned grids with cut cells at 
the wall. We are using a new wall model recently developed by Steven Allmaras. The derivation of the 
model is described in Appendix C of this paper, written by Allmaras. We repeat the final equation here for 
convenience, 

u + {y + ) = B + cilog((y + +ai) 2 + b\) - c 2 \og((y + + a 2 ) 2 + bl 
— C3 arctan(j/ + + cq, 61) — C4 arctan(j/ + + a2, 62), 


where the constants are given in the Appendix. Several other researchers 13, 36 have used Spalding’s composite 
formula, an algebraic relation bridging the viscous sublayer and the log layer within a fully developed 
turbulent boundary layer, as a wall function. For completeness we give Spalding’s formula here too, 


y + = u + 


-kB 


e ku+ - 


1 — ku + — t( ku + ) 2 — i(fcw + ) 3 ^ 


( 13 ) 



Figure 5. Comparison of new SA wall model with 
Spalding’s formula. Most of the difference is in the 
transition region. 


Equation ( 12 ) (henceforth called the SA wall 
model) has the advantage over Spalding’s formula of 
matching the profiles actually computed in the flow 
field by the Spalart- Allmaras turbulence model. This 
is seen in the computational examples of Section C, 
where the computed profiles match the SA model 
through the transition region. Figure 5 shows the 
transition region, where the two formulas differ the 
most. In this figure, Spalding’s formula was evaluated 
using B = 5.033 in ( 13 ) to match the asymptotic be- 
havior of the SA model, which requires this value for 
the constant shift. Spalding’s formula simply bridges 
this buffer zone with a functional form. Furthermore, 
(12) is more computationally efficient, since it is a for- 
ward evaluation for u + given ?/ + , rather than needing 
a nonlinear iteration as ( 13 ) does. (However Kalitzin 
et al. 15 have a nice approach that avoids iteration us- 
ing pre-computed tables of inverses) . By using a wall 
model that actually matches the background flow we 
can avoid the kinetic energy mismatch, described e.g. 
in Sondak and Pletcher. 37 


2 . Coupling 

Since the model is given in wall coordinates (t/ + , ?x + ), y + = u + = an estimate of the friction velocity 
u T (a.k.a. du/dy\ wa u) is needed for each cut cell. Wall functions have a well-known sensitivity to the location 
of the first point away from the wall, 34 and this issue must be directly addressed for the wildly irregular cut 
cells of an embedded-boundary mesh. Following Capizzano 18 we regularize the distance by using points a 
fixed distance h away from the wall instead of the using the cut cell centroid. 

We fit the one-dimensional model in (12) using the boundary condition of zero velocity on the wall and 
the solution at a fixed point F (to use the same terminology 18 ) located a distance h from the boundary in 
the normal direction. Let the point F be in cell D; the approximate solution there is obtained using linear 
reconstruction from D’s centroid. The velocity is then rotated into qtang , the tangential velocity, using the 
directions defined by the portion of the boundary in cut cell C. A Newton iteration is used to find the friction 
velocity u T that transforms the point (h, qtang{F )) into (h + , qta ng (F)) lying on the SA wall model curve (12). 
Gradients of the model provide the viscous flux needed at the wall for the finite volume scheme. The model 
also gives values of the tangential velocity needed to compute the difference at the cut faces (after rotating 
back to Cartesian velocity components u and v). Finally, as with the quadratic we replace the cut-cell least 
squares gradient with a model gradient (and set the tangential component to zero). This coupling strategy 
makes use of the fact that the friction velocity is constant through the inner portion of the boundary layer. 
The procedure is illustrated in Fig. 6. 
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Figure 6. Illustration of construction for wall model coupling via constant height “forcing” points through the 
centroids of the cut cells. 


Currently we take the distance h — 1.5 x min (dx,dy). Since the cell sizes in the coordinate directions 
are approximately equal, this distance is slightly larger than a diagonal in a two dimensional cell, and so 
is sufficient to situate the point F out of the cut cell and into the neighboring cell. This procedure does 
impose some restrictions on the mesh resolution, since the point F must lie in the logarithmic region of the 
boundary layer to get valid estimates of the friction velocity. Eventually this can be used as a criteria for 
adaptation, or alternatively the location of the point F can vary depending on the local flow conditions. 

Section IV. C shows results of coordinate- aligned turbulent flat plate simulations with Reynolds number 
ReL = 5 x 10 5 using both approaches - integrating down to the wall and using the wall model. The latter 
needs substantially less resolution. We will also show a comparison to the more difficult case where the flat 
plate is rotated with respect to the mesh. The wall model computation needs an extra level of refinement 
in the non-coordinate-aligned case. These preliminary studies show savings commensurate with the savings 
found in body- fitted grid implementations of wall models. 16,38 


III. Multigrid Acceleration and Time-Stepping 

Our steady-state solution strategy uses an FAS multigrid method with a multi-stage Runge-Kutta 
smoother applied on each multigrid level. It is common for viscous flows on highly anisotropic grids to 
cause convergence problems for multigrid in getting to steady-state. However, since the Cartesian cells in 
this work are essentially isotropic, appropriately strengthened multigrid operators should yield Navier-Stokes 
performance on a par with our inviscid multigrid algorithm. 22 In RANS simulations, the equation for the 
turbulence model also needs time advancement until a steady state is reached. Since this equation has 
stiff source terms that may restrain the rest of the system, we use a separate timestep for this equation 
as described below. The system is solved fully coupled, with the Runge-Kutta smoother operating on all 
equations, but using one timestep for the flow equations and another for the turbulence model. 


A. Time-Stepping 

We compute directional timesteps for the flow equations following the usual advect ion- diffusion model. 


A t x = 
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where c is the local sound speed and 
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The timestep within the cell is a harmonic blend of the directional AFs. 
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Following Pierce, 39 the use of a separate timestep for the turbulence model in a fully coupled time- 
advance scheme amounts to a diagonal timestep formulation, essentially applying a scalar preconditioning to 
the turbulence model. The timestep for the SA working variable, may be written concisely by replacing 
K v with K v SA 

^ = 2(v + ,)(l^) (17) 

resulting in directional timesteps 


A t x 


fox 


M + 1 + 


K.SA > 
hx 


Aty — 


hy 

H + f + 


Ku SA ' 

ky 


(18) 


The factor | is an ad-hoc inclusion for robustness that prevents the hyperbolic contribution from vanishing. 39 
The directional timesteps for the turbulence model are again blended via the harmonic mean in (16). 

The source term (“Production - Destruction”) is treated implicitly by penalizing the timestep when its 
Jacobian, is negative. 30,39 


Atimr, — 


At 

1 — At min(0, Q„) 


(19) 


B. Second-Order Coarse Grid Operator 

Inviscid multigrid solvers commonly use a first-order spatial differencing on coarser grids. This saves the 
expense of computing gradients, since the coarse grid solution does not affect the final solution. Without 
gradients, the only geometric information required is the area of the cut faces, the cell volumes, and the wall 
normal vector in each coarse cell. By contrast, a second-order coarse operator needs solution gradients to 
perform linear reconstruction to coarse grid cell faces. Near the wall, this requires the cut-cell centroids, cut- 
face centroids, and surface (wall) centroids on the coarser grids. Most of this information is easily computed 
from the fine grid and can be done concurrently with coarse mesh generation. For example coarse grid cell 
centroids are weighted averages of all the fine grid cells that restrict to that cell. 

The most difficult part of the coarse grid generation is organizing split cells on-the-fly during mesh 
coarsening. Split cells are cut Cartesian cells that may be split into multiple control volumes by the geometry. 
Figure 7 shows two interesting examples. Cut cells on the fine grid may agglomerate into split cells on the 
coarse mesh with multiple control volumes (Fig. 7a). Alternatively, split cells on the fine grid can yield a cut 
(but unsplit) coarse cell (Fig. 7b). Specific situations like these may be quite complex in three dimensions 
since they can involve any number of control volumes. Nevertheless, all such cases can be resolved using a 
simple integer matching algorithm. This algorithm scans for common triangles on sorted lists of intersected 
triangles maintained by each cut cell, along with the face lists that connect cut cells with their neighbors. 
Since it is based on integer comparisons of sorted lists, this matching procedure is very fast and robust. 



Figure 7. (a) Four fine cut cells that make one coarse split cell containing two control volumes, (b) A fine 

split cell, a cut cell, and two full cells that make up one coarse cut (but not split) coarse cell. 


C. Linear Prolongation Operator 

The other major change in the multigrid algorithm is the formulation of the prolongation operation. Multigrid 
theory shows that higher-order derivatives need better prolongation than the piecewise constant approach 
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typical of inviscid flow with cell-centered schemes. 40 Our restriction operator, a volume-weighted average of 
the solution and residual on the fine grid, is already second-order, and a second order coarse grid discretization 
was also implemented. With gradients now defined on the coarse levels it is easy to use a linear prolongation 
to bring the change in the solution back to the finer levels. Two copies of the solution were already needed - 
the initial restriction to the coarse grid and the modified solution after recursively smoothing on the coarser 
levels. Thus, in addition to the computational expense, there is some minor additional memory overhead 
from storing gradients on the coarser levels. 

Figure 8 compares convergence of the multigrid 
scheme before and after improvement for the Navier- 
Stokes equations (1). The black curve shows the orig- 
inal formulation without coarse gradients and using 
only piecewise constant prolongation. The blue curve 
shows the improved operator using a second-order 
coarse grid scheme combined with a linear prolonga- 
tion operator. While this particular data is for the 
laminar (Rcl = 5000) flow over a NACA 0012 air- 
foil at Mqo = 0.5 (presented in Section IV. B) the 
improvement is characteristic of the scheme’s per- 
formance. In this example, six levels of multigrid 
were used with a W-cycle consisting of one pre- and 
one post-sweep of a 5 stage Runge-Kutta smoother 
on each level. With linear prolongation and coarse 
grid gradients, the scheme remained stable with a 
CFL number of about 1.4. In contrast, the unmodi- 
fied scheme in Fig. 8 had a maximum stable Courant 
number of only about 0.1. Note that the improved 
convergence and stability required the combination 
of both linear prolongation and the improved coarse 
grid spatial scheme - using either alone was insufficient to improve multigrid performance. Figure 8 shows 
that when used in combination we can achieve convergence behavior nearly on a par with that of the base 
inviscid scheme previously reported. 22 

Finally, Fig. 9 shows multigrid convergence be- 
havior of a fully turbulent RANS example with the 
Spalart-Allmaras turbulence model. The fully cou- 
pled time advance is using the diagonal timestep de- 
scribed earlier in combination with the second-order 
coarse discretization and linear prolongation opera- 
tor. The case shown corresponds to the M ^ = 0.676 
fully turbulent flow over an RAE 2822 airfoil at a 
Reynolds number of 5.7 x 10 6 presented in Section 
IV. D (medium mesh, AGARD Case 1) computed 
with 5 levels of multigrid. Cells at the wall are 
about three times finer than in the NACA exam- 
ple of Fig. 8. While the additional stiffness intro- 
duced by the turbulence model has clearly impacted 
convergence, the scheme is still delivering roughly 9 
orders-of-magnitude reduction in the L\ norm of mo- 
mentum in under 1000 (fine-grid) multigrid cycles 

for both the flow and turbulence model. This cor- 

i, u . . i , r , Figure 9. Multigrid convergence for fully turbulent 

responds to a multigrid convergence rate of about „ „ J n ___ , 

1 ... flow over an RAE 2822 airfoil at = 0.676 and 

0.98. While this is noticeably slower than the rate of Reynolds number 5.7 X 10 6 (see Section IV. D) using 
roughly 0.92 shown in laminar case (Fig. 8), it is still a five level multigrid W-cycle. 

fast enough to support our numerical investigations. With approximately 115,000 cells in the mesh, this 
example took about 5 minutes on a current generation desktop computer. The literature on block- Jacobi 
preconditioning and matrix time- stepping offers many suggestions for further improvement of these results. 


Figure 8. Multigrid convergence for NACA 0012 with 
Moo = 0.5, and Re^ = 5000 at zero angle of attack. 
The figure compares the use of linear prolongation and 
coarse grid gradients with piecewise constant prolon- 
gation and no coarse gradients. Six levels of multigrid 
were used. 


FMG 
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IV. Computational Examples 


In this early investigation, computational results are principally confined to verification and validation 
exercises. Our goal is to investigate the accuracy of the viscous discretization operator on cut-cell Cartesian 
meshes and examine fundamental issues surrounding stencil irregularity, wall modeling and resolution re- 
quirements. We do not focus on performance, since we have not yet seriously begun to attack the cell count 
or meshing efficiency issues. 


A. Laminar Flat Plate Boundary Layer 

The first test is a two-dimensional flat plate boundary layer for M ^ = 0.5 and ReL = 5000. As in Section 
II the plate is oriented 15° to the mesh. The plate starts at x = 0 in a domain with x spanning -9 to 20 and 
y from -3 to 30. The HLLC Riemann solver 41 was used, with wave speeds from Batten et al . A 2 As discussed 
in Section II. 2, these laminar computations all used quadratic reconstruction in the wall- normal direction to 
obtain the necessary derivatives in the cut cells. 

Figure 10 shows a sketch of the problem setup (top left) along with profiles of both tangential and normal 
velocity which were extracted from solutions with and without limiting. The velocity profiles are taken at 
3 stations corresponding to Re x = 5000, 10000, and 50000 along the plate, plotted in similarity coordinates. 
As expected, the tangential velocity profiles collapse on each other in these variables. The normal velocity 
profiles show some effects of the plate leading edge and the mesh resolution. The calculation on the left did 
not use limiters, and shows some viscous overshoot at the first two stations. Velocity profiles at the right of 
Fig. 10 used the van Leer limiters and have no viscous overshoot. Isotropic Cartesian cells were used with a 
cell size at the wall of h = 5.9 x 10 -3 , giving approximately 13, 17, and 40 cells in the boundary layer at the 
three stations respectively. These resolution requirements are similar to standard second-order finite volume 
codes. 43 Note that since the profiles are taken in the direction normal to the wall and not a coordinate 
direction, each profile intersects both x and y grid lines along the way, so the number of symbols does not 
correspond to the number of cells. 


No Limiter 




van Leer Limiter 



Figure 10. Three profiles at Re x = 5 K, 10 K and 50 K from the flat plate boundary layer for Moo = 0.5, ReL = 5000. 
The plate is rotated 15° to the mesh, shown on the left. In the middle figures, the profiles were not limited; 
on the right the van Leer limiter was used. 
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Figure 11. Skin friction for flat plate at 15° to the mesh using the quadratic is well predicted by Re x = 100. 

Figure 11 shows the skin friction distribution for this case, discussed in detail in Section II. Cf is well 
predicted by Re x = 100, and compares nicely with coordinate-aligned results. 43 Note the smooth Cf 
distribution provided by the quadratic wall-normal reconstruction. 

B. NACA 0012 

Another test is laminar flow about a NACA 0012, also at M ^ =0.5, Reynolds number 5000, and zero angle 
of attack. Unlike the flat plate, in this case the stencil of the quadratic used in each cut cell is not in a 
fixed direction, but rotates around the airfoil. The results are stable regardless of the direction of the stencil. 
Figure 12 shows an overview of the discrete solution along with some examples of the quadratic’s stencil near 
the leading and trailing edges. What sometimes appears to be 3 connected points in the stencil are really 
2 stencils of 2 points each for adjacent cut cells. Figure 12b shows the velocity magnitude, computed on a 
grid with leading edge resolution of 0.0016c obtained with 10 levels of refinement along the airfoil. Figures 
13a and b show the pressure and skin friction coefficients around the body. Again the skin friction is smooth 
along the airfoil. Separation here occurs at x = 0.816, very close to the mesh converged values found in the 
literature. 43 


Leading edge 




Trailing edge 

(a) 


NACA 0012 



Moo = 0.5 
Rei_ - 5000 


(b) 


Figure 12. (a) Cells used in quadratic stencil for computing wall fluxes for NACA0012 example, (b) Velocity 
magnitude. 
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-U.D 


0.15 



x/c ' ’ x/c 

Figure 13. Cp (left) and Cf (right) for the results in Fig. 12b for flow around a NACA 0012. 


C. Turbulent Flat Plate 

While turbulent boundary layers are thicker than laminar boundary layers, the higher wall shear stress creates 
a greater demand for resolution near the wall. In this example we use the RANS equations to simulate a flat 
plate with free stream Mach number 0.5 and fully turbulent flow at ReL = 5 x 10 5 . We will first compare 
wall treatments with and without the wall model for a plate that is coordinate-aligned. We then compare 
the use of the wall model in the aligned case with that of a plate angled at 15° to the mesh. All cases uses 
the Spalart-Allmaras turbulence model; the cases with the wall model use the SA wall model. 

The baseline case of integrating the governing equations down to the wall using 13 levels is shown in 
Fig. 14. The tangential profiles are taken from 4 stations along the plate. They are all well-resolved and 
show no viscous overshoot. We also compare the velocity profile at x = 10 corresponding to Re x = 5 x 10 6 
in plus coordinates with the SA wall model. With an initial cell at a y + = 8.8 the results match well, and 
are virtually identical to 14 level results (not shown). Skin friction is shown on the right, compared to the 
asymptotic formula for a flat plate with zero pressure gradient taken from White. 44 The solution is not noisy, 
since the wall is aligned with the mesh and all cell centroids are the same distance from the wall. After 13 
levels of refinement at the wall, the Cartesian cells have a mesh spacing of h = 3.6 x 10 -4 . 

The computational results using the SA wall model are shown in Fig. 15. The figure compares wall 
model results in both the aligned and non-aligned case. With the wall model, accurate results are achieved 
on an aligned mesh with only 10 levels of refinement. No attempt has yet been made to refine the mesh in 
a solution-adaptive manner; it is only refined next to the flat plate. Nevertheless we can report that our 




Figure 14. Coordinate-aligned turbulent flat plate with ReL = 5 X 10 5 , Moo = 0.5. Left figure shows well resolved 
profiles with no viscous overshoot. Calculations were done integrating down to the wall using the quadratic 
on a 13 level mesh. Middle figure shows the velocity profile at Re x = 5 X 10 6 plotted in wall coordinates and 
compared to the SA formula for a zero pressure gradient flat plate. For this profile the first cut-cell centroid 
corresponds to a y+ « 8.8. Right figure shows skin friction compared with asymptotic formula from White. 44 
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non-optimized 13 level mesh had 1.6M cells, and the 10 level mesh had 180K cells. 

The computation is much more difficult on the 
non-aligned plate, since the difference scheme for 
the cut cells degenerates to first-order. Results on 
the flat plate oriented at 15° to the mesh needed 
one finer level of mesh refinement. In addition 
to the two results with the wall model, Fig. 15 
also includes the results of integrating down to 
the wall. The aligned 13 level result integrat- 
ing to the wall (y + = 8.8), the aligned 10 level 
wall model calculation (y + ^ 70), and the non- 
aligned 11 level results (y + « 35). all agree quite 
well for small y + . In the wake region the over- 
shoots in the solution clearly shows the effect of 
the larger truncation error when the streamwise 
direction is not aligned with a coordinate direc- 
tion. Additional mesh refinement would certainly 
help in this region of high curvature around the 
knee of the profile. Alternatively, future research 
can try to improve the accuracy of the discretiza- 
tion to fourth-order, an especially attractive idea 
on Cartesian meshes. (Note that the first point in the plot of y + versus u + is taken from the forcing point 
F, and so does not necessarily correspond to the finest mesh spacing reported at the wall.) 

D. RAE Airfoil 

To assess performance and resolution requirements 
of the numerical scheme and wall model for high- 
Reynolds number turbulent flows, we consider three 
cases using the RAE 2822 airfoil. Mesh conver- 
gence studies were conducted on a sequence of three 
meshes for AGARD Cases 1, 6 and 10. 45 The coarse 
mesh in this sequence has a wall spacing of 0.1% of 
the airfoil chord, c, and a total of about 50k cells. 

For these cases, this wall spacing gives values of y+ 
over the airfoil largely in the range 200-350. The 
medium mesh has a wall spacing of 0.05%c and a 
total of ~115k cells with y + of 100-170. The fine 
mesh has a wall spacing of 0.025%c and ~150k cells 
with y+ values largely between 50 and 80. The 
medium and fine grids were constructed by sub- 
division of the near-wall cells of the coarse mesh, 
thus the three grids are nesting, and the domain 
extended a distance 30c from the airfoil. No cir- 
culation correction was applied in the simulations Figure 16. Mach contours of flow around RAE 2822 
and all simulations were run fully turbulent. airfoil at AGARD Case 1 conditions. = 0.676, Re L = 

5.7 x 10 6 and Cl = 0.566 ( a = 1.85°) using the SA wall 
model. 

1. AGARD Case 1 

The first case we consider is AGARD Case l. 45 This is a subcritical flow at = 0.676 at a Reynolds 
number of Re c = 5.7 x 10 6 . Experimentally, the target lift coefficient of Cl = 0.566 was achieved at an 
uncorrected a = 2.4°, however most simulations match this value with about half a degree less incidence 
angle. 46 

Figure 16 illustrates this flow through Mach contours on the fine mesh. The simulation matched the 
experimental value of lift at a = 1.85° and subsequent simulations on the coarse and medium meshes 
were run at the same incidence angle. Figure 17 displays the surface pressure coefficient and skin friction 




Figure 15. Non-aligned turbulent flat plate using the wall 
model on an 11 level mesh, compared with the aligned 
cases, both using the wall model and using the quadratic. 
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Figure 17. Surface pressure coefficient and skin friction on RAE 2822 airfoil for AGARD Case 1 conditions 
using the SA wall model on coarse, medium and fine meshes with experimental data from Cook et a/. 45 


distribution for all three meshes and contains the experimental data for comparison. Agreement between 
simulation and experiment is on a par with other published solutions, 46 and the presence of simulation data 
from three meshes makes it possible to evaluate the level of mesh convergence. The pressure distribution (left 
of Fig. 17) shows good capturing of the suction peak, the rooftop region and the characteristic “duck-tail” 
in the highly-loaded aft-portion of this airfoil. Overall, the pressures are remarkably invariant with mesh 
resolution, indicating that these pressures are essentially mesh converged even on the coarsest grid. Since 
the aft recovery region depends upon prediction of the displacement thickness, it is worth noting that even 
the coarsest mesh seems to accurately capture this feature despite local y + values of approximately 200 on 
this grid. 

Since it is a derivative quantity, prediction of skin friction is more challenging (Fig. 17 right). The 
experimental drag value was measured at 85 counts. The simulations are a bit higher, showing Cd = 
{0.0099, 0.0097, 0.0099} on the coarse, medium and fine meshes respectively. Despite similarities in the skin 
friction distributions, and the consistency in the integrated drag, there are significant differences between the 
profiles. At the trailing edge, where the boundary layer is the thickest, local values on the three meshes are 
in good agreement with both each other and experimental data. Following upstream from the trailing edge, 
differences become apparent. As the boundary layer thins, the coarse mesh results start to diverge. The 
fine and medium meshes agree over most of the surface, however upstream of about 30%c, even the medium 
mesh appears to be insufficient as it starts to diverge from both the fine mesh solution and the experimental 
data. Near the leading edge, only the fine mesh remains credible and shows good agreement with the data 
upstream of 20% chord. As in the experimental dataset, skin friction data in this figure was normalized 
using the dynamic pressure at the outer-edge of the boundary layer to facilitate direct comparison. 


2. AGARD Cases 6 & 10 

The final two examples both consider supercritical flow over the same RAE 2822 airfoil section and are 
characterized by increasingly strong shock-boundary layer interactions. AGARD Case 6 was tested at = 
0.725, Re c = 6.5 x 10 6 and produced a measured lift coefficient of Cl = 0.743 with an uncorrected incidence 
angle of a = 2.92°. Case 10 was tested at the same lift coefficient, but under somewhat more aggressive 
conditions: = 0.750, Re c = 6.2 x 10 6 . Cl = 0.743 was measured at an uncorrected incidence of 

a = 3. 19°. 45 Over the past 30 years, these two cases have been the subject of countless numerical experiments. 
As in Case 1, most numerical simulations achieve the experimental values of lift in these cases with around 
half a degree less than the experimentally reported incidence angle. 46 

Figure 18 illustrates these two flows through isoclines of Mach number with case 6 at the left and case 
10 at the right. The fine mesh simulation of Case 6 matched the experimental lift coefficient of Cl = 0.743 
at a = 2.44° while the simulation in Case 10 matched this same value at a = 2.54°. As before, simulations 
on all meshes were performed using these values. Despite identical lift coefficients, the higher Mach number 
in the Case 10 flow produces a stronger shock which is shifted noticeably aft. 

Figure 19 contains a comparison of surface pressure distributions for these two transonic cases. Agreement 
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Figure 18. Mach contours of flow around RAE 2822 airfoil at AGARD Case 6 (left) and Case 10 (right) 
conditions using the SA wall model on the fine mesh. Case 6: = 0.725, Re c = 6.5 x 10 6 @ Cl = 0.743 

(■ a = 2.44°). Case 10: = 0.750, Re c = 6.2 x 10 6 @ C L = 0.743 (a = 2.54°). 


with the experimental data is comparable to others in the literature for these cases. 30,46,47 In Case 6, the 
computed pressures on the lee side do not dip as much following the suction peak as the experimental data, 
and consequently the shock is slightly forward to yield the same lift coefficient. In Case 10, the numerical 
solution places the shock slightly farther aft than the experimental and shows small discrepancies in the 
aft loading. Similar differences have been noted by numerous authors. 30,46,47 As in the earlier subcritical 
example, both of these cases show little variation of the pressure distribution as the mesh is refined, and even 
the coarse mesh with a wall spacing of 0.1%c and only 50k cells seems sufficient to produce mesh- converged 
pressures. The shock position in both of these examples is critically dependent on the displacement effect 
of the boundary layer and the shock-boundary layer interaction and the performance of the coarse mesh is 
noteworthy in this regard. 

Figure 20 displays the the evolution of Cf for both Case 6 and Case 10 as the mesh is refined. The profiles 
display similar behavior to that discussed in the subcritical example earlier (Fig. 17). Integrated drag on 
the fine mesh was 133 counts which compares well with the experimental value of 127 counts. 45 As before, 
agreement is strongest where the boundary layer is the thickest. Tracing upstream from the trailing edge 




Figure 19. Surface pressure coefficients on RAE 2822 airfoil at AGARD Case 6 (left) and Case 10 (right) 
conditions using the SA wall model on coarse, medium and fine meshes with comparison to experimental data 
from Cook et al . 45 . 
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Figure 20. Skin friction distributions on RAE 2822 airfoil at AGARD Case 6 (left) and Case 10 (right) 
conditions using the SA wall model on coarse, medium and fine meshes with comparison to experimental data 
from Cook et a/. 45 . 

reveals that the coarse mesh starts to diverge at about 80%c. The medium mesh shows discrepancy upstream 
of 30%c, and only the finest mesh accurately predicts the experimental data near the leading edge where the 
boundary layer is the thinnest. This case is characterized by the strong shock-boundary layer interaction 
on the suction side. While results in the literature are mixed, many authors report a small recirculation 
bubble following this shock. 46 The present simulations are showing very small, but still positive, values of 
skin friction indicating only incipient separation. Detailed investigation of the analytic wall model’s behavior 
near the shock-boundary layer interaction remains an important topic of investigation. 

V. Conclusions 

In this paper we outlined initial development of a method for solution of the Reynolds Averaged Navier- 
Stokes equations on embedded-boundary Cartesian meshes using a wall model in the cut cells. First we 
developed accurate discretizations based on recentering for use in the highly irregular cut cells. We also 
introduced an analytic wall model for turbulent flow simulations based on a limiting solution of the Spalart- 
Allmaras turbulence model which provides accurate representation of attached flows in the laminar sub- 
layer, buffer and log-layer regions of the boundary layer. This model has distinct advantages over standard 
Spalding-based wall functions since it automatically matches the turbulence model used away from the wall. 
Additionally, it can be written as a forward evaluation for the velocities as a function of distance and does 
not require iterative solution. We use an enhanced multigrid algorithm that includes geometry and gradients 
on coarser grids to achieve multigrid convergence nearly as good as in the inviscid case. Finally, we showed 
that the resulting scheme can provide accurate pressure distributions for high Reynolds number flow over 
subsonic and transonic airfoils with wall spacings of about 0.1 %c - even when the pressure distribution is 
strongly dependent on the displacement of the boundary layer. Drag values for attached flows are reasonable 
at this same resolution, but detailed skin friction distributions require at least one or two additional cell 
refinements. 

Resolution requirements of the current method are driven by the assumptions underlying the analytic wall 
model. The new SA wall model, like a Spalding-based wall function, is still based on a simple diffusion model 
of the near-wall flow and so it cannot accurately predict values outside the log layer. Both the non-aligned 
turbulent flat plate and 2D airfoil examples showed that accurate local skin friction values required resolution 
to the log layer. This is currently the limiting factor in determining how coarse a Cartesian mesh can be 
used. This observation provides strong motivation for work towards a more comprehensive PDE-based wall 
model which retains more physics from the governing equations. 

The present results are an important first step, but clearly a long list of questions remain. Our immediate 
goal is development of a more complete, non-analytic, wall model which includes streamwise momentum and 
pressure gradient terms that become important outside the log layer and in separated flows. Robust conver- 
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gence and coupling of such a model with the underlying Cartesian method is likely to be another challenge. 
A successful three-dimensional implementation will clearly depend on the wall model to mitigate the resolu- 
tion requirements that would lead to untenable cell-counts. Nevertheless, the simplicity and automation of 
Cartesian mesh generation continue to provide persuasive motivation for further algorithmic development. 


Appendix A: Numerical Discretizations of Second Derivatives at Cut Cells 


The development of an accurate discretization for 
the Navier-Stokes equations in a finite volume context is 
much more delicate than for the inviscid equations. Two 
new stencils need to be developed to compute second- 
derivative terms at embedded boundaries - one at the 
cut faces between cut cells, and the other to compute the 
cut-cell boundary conditions at a solid wall. We consider 
three possible discretizations for an elliptic equation to 
evaluate their accuracy on a problem with an analytic 
solution. They were chosen because they are nominally 
second-order, conservative and can be implemented in 
a finite volume framework. These tests do not rely on 
special properties of elliptic equations. 

Consider the model Poisson equation in two dimen- 
sions 

v 2 0 = 



Figure 21. Domain for elliptic model problem in 
Table 1. 


/ 


( 20 ) 


with exact solution (ft(x,y) = sin(x) exp(?/). The domain is a quadrilateral with three non-coordinate-aligned 
sides, as shown in Fig. 21. Three of the four sides use Neumann boundary conditions, with Dirichlet on the 
fourth. This is in contrast to the Euler equations, where only the value of the solution itself is needed at the 
cell edge. At interior cells, the standard central difference approximation for the first derivative, resulting in 
the second-order accurate 5-point Laplacian is used. Three alternatives are tested for the cut cells. 


Recentering: The first method uses a “recentering” idea 48,49 illustrated in Fig. 22a. Recentering 
operates on the integral cell average, which for a second-order finite volume scheme can also be used for 
the pointwise value of the state vector at the cell centroid. Within each cut cell, a least squares gradient is 
reconstructed using the solution from each cell’s stencil of face neighbors. The gradient is used to reconstruct 
the solution from the cut-cell centroid to a line normal to the face through the face centroid (see Fig. 22a). 
This is done on both sides of a cut face, so that the derivative can be approximated by a simple difference 
in the face normal direction. Using the notation of Fig. 22a, with cell centroids at A and (7, we recenter to 
locations B and D , defining 


(ft X 


ft>B — ftp 
Xb — %D 

((ft A + V0A • (Iba) - {(ftc + V0C • dpc) 
X B - X D 


( 21 ) 


where cIba is the vector from A to B, and similarly for doc- Note that this divided difference is nominally 
first order accurate, since it is not centered about the face centroid. Changing the recentering step to use 
locations that are equally spaced about the face centroid, so that a centered difference approximation can 
be used, is no more accurate than the one-sided scheme without using a more accurate reconstruction, so we 
do not include those results here. 

Johansen-Colella: The second discretization uses the Johansen-Colella framework. 50 In this approach 
the solution <ft is thought of as being located in the center of the original uncut cell, not at the cut-cell 
centroid. First derivatives are then easily computed at the midpoint of an uncut face. The flux however 
is needed at the cut-face centroid, illustrated in Fig. 22b at the point P. This is easily obtained by linear 
interpolation along the edge. 
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(a) Recentering: The solution is 

reconstructed from cell centroids (A 
and C) to new points in each cell (B 
and D) on the line perpendicular to 
the midpoint of the cut-cell edge. A 
simple difference (B-D)/Ax approxi- 
mates the derivative. 



(b) Johansen-Colella: The solu- 
tion is considered to be located in 
the center of the uncut cell. The flux 
at the centroid of the cut face, P, is 
linearly interpolated from the uncut 
fluxes (A21 and A43 ) . 



(c) Polynomial Fit: A least 

squares polynomial that is linear in 
x and quadratic in y is fit using 6 
neighboring values. The polynomial 
is differentiated and evaluated at the 
point P to approximate the flux. 


Figure 22. Three methods for computing the derivative at a cut-cell edge. 


Polynomial Fit: The third discretization is similar to a method in Ye et al. 51 At each face, a least- 
squares polynomial that is linear in the face normal direction (e.g. the x direction for an x face) and quadratic 
in the transverse direction (e.g. y direction for an x face) can be constructed using 6 neighboring centroidal 
values. We use the solution to the left and right of the face, the boundary values in those cells, and the 
values one cell further away, illustrated in Fig. 22c. The polynomial is differentiated and evaluated at the 
cut-face centroid to obtain the flux. 


In all three methods the wall flux was computed using the same discretization (modulo the location of 
the variables). In finite volume form the flux at the wall (d(p/dn) is needed. For the sides with Neumann 
boundary conditions this becomes an evaluation of the boundary condition. On the Dirichlet side we use 
a one-sided discretization, taking the cell average, subtracting the Dirichlet boundary condition evaluated 
at the point normal to the centroid (or cell center, in the case of Johansen-Colella), and dividing by the 
distance between them. 


Recentering Johanson-Colella Polynomial fit 


Size 

1 norm 

max norm 

1 norm 

max norm 

1 norm 

maxnorm 

9x9 

38.53 

1.47 

24.64 

1.39 

46.93 

4.76 

18 x 18 

10.10 

.40 

7.18 

.48 

10.66 

1.44 


conv. rate 

3.8 

3.7 

3.4 

2.9 

4.4 

3.3 

36 x 36 

2.60 

.11 

1.81 

.12 

2.64 

.60 

conv. rate 

3.9 

3.6 

4.0 

4.0 

4.0 

2.4 

72 x 72 

.66 

.03 

.45 

.03 

.65 

.13 

conv. rate 

3.9 

3.7 

4.02 

4.0 

4.0 

4.6 


Table 1. Results of solving the Poisson equation V 2 </> = / on an irregular domain using 3 different discretizations 
for the cut cells. Data corresponds to Fig. 23. 


Figure 23 and Table 1 report the error in computing the solution to (20). Results include both the 
maximum norm of the error, HeH^ = max^- 1 4 > exacts ~ computed . . |, and the LI norm of the error, ||e||i = 
JA • Aij | (j) eX act i j — 4> computed XJ |> where the Aij are the cell areas. (Note that the LI error is not normalized by 
the domain size, so it is larger than the maximum norm errors.) All methods show second-order convergence 
in the LI norm. The convergence is not entirely smooth, since cut cells do not have a smooth asymptotic 
expansion for the error. 
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We observe that the Johansen-Colella scheme has 
slightly better L\ performance, and that the polyno- 
mial fit is somewhat less accurate (especially in the more 
finicky maximum norm). This conclusion is representa- 
tive of a variety of test cases we ran. Since the accuracy 
of the recentering approach and Johansen-Colella is sim- 
ilar, and the former fits better into our existing finite vol- 
ume cut-cell framework, the recentering approach was 
chosen for implementation of the Navier-Stokes equa- 
tions. These results come with a caveat, and need to 
be interpreted carefully. There are terms in the Navier- 
Stokes equations that are not included in this model 
problem. In practical settings, these remaining terms 
can dominate the error, and lead to non- convergence if 
not handled carefully. In addition, we are not looking 
at positivity of the stencil, which has been the focus of 
several other studies. 2,26 Although positivity ensures a 
maximum principle, blindly insisting on it may rule out 
potentially useful discretizations. 3 



Figure 23. Convergence of 3 methods for comput- 
ing cut-cell fluxes for the model Poisson problem 
V 2 b = /, corresponding to data in Table 1. 


Appendix B: Stability of Quadratic Boundary Condition 

This analysis examines the effect on the allowable time step of using a quadratic instead of a linear 
interpolant to discretize the boundary condition for a model problem. We will use GKS theory 29 and 
consider the heat equation u t = u xx on [0,1], with u(0) = u( 1) = 0. Using forward Euler in time and central 
differencing in space gives = u™ + 1 —Uj)/h — ( Uj — Uj-i)/h), where the flux term is written in 

finite volume form. The boundary condition u(x = 0) = 0 is usually implemented as .b(u\ + uq) = 0 where 
uq is a ghost cell. A simple calculation shows that this boundary conditions does not reduce the stability 
limit A = At/ti 2 = .5 of the initial value problem. 

Using divided difference notation and the Newton form 52 of the polynomial, the quadratic interpolant of 
the three values u b ,ui and u 2 can be written 

p 2 (x) = [u b \ + [Ub,ui\x+ [u b ,ui,u 2 \x (x-xi) 

where we have written u b to explicitly represent the boundary value at x = 0 even though it is zero. The 
left flux at the first cell will be discretized using p' 2 at the left cell edge, P2 f (0) = {—u 2 + 9ui)/3h so that the 
equation for the update of u\ is 


up 1 = u™ + — ({u 2 - Ui)/h -p' 2 ) (22) 

It is sufficient to consider the stability of the left half plane problem, and look for solutions of the form 
u'j = k? z n . Substituting this into the difference schemes gives the characteristic equation 

2 = 1 + A(/c — 2 + l/«) (23) 

where we are interested in solutions with k < 1. The characteristic equation for the boundary condition (22) 
is 

2 = 1 + ^(4k-12). (24) 

o 

Equating (23) and (24) gives an l 2 solution k = 3 — 2a/3- 

Substituting this root into (24) and solving for z < 1 gives A < a/ 3/4 ~ .43. This is only a small 
reduction in the stability limit for this model heat equation. Since viscous terms in the Navier-Stokes 

3 For example, the 4th-order Laplacian derived using Richard extrapolation by combining the 3-point h-sized stencil and the 
2 h stencil has negative coefficients on the points 2 h away from the center yet is symmetric positive definite. 
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equations only reduce the allowable time step by a small fraction, the impact on the time step in practice 
would be correspondingly less. 


Appendix C: Analytic Law of the Wall Solution 
for the Spalart-Allmaras Turbulence Model 

by Steven R. Allmaras 


We make the usual assumptions for law of the wall analysis: incompressible, zero pressure gradient, 
constant outer (edge) velocity, advection terms are negligible and d/dx d/dy, where x is streamwise and 
y normal to the wall. With these assumptions, u = u(y) and u = u(y), and the r-momentum equation 
reduces a statement of constant total shear stress, 

du 

— == const = u 
dy 


(25) 


d 

dy 


. . du 

(v + vt) ~r 

dy 


= 0, 


iy + vt) 


where u T is the wall shear stress velocity. With these same assumptions, Spalart-Allmaras (SA) reduces to, 


1 d 

a dy 


, dv 


Q>2 

a 


du 

dy 


1 2 _ 

+ c b i (1 - f t2 ) Sv - (c wl f w - ^f/ t2 ) 


r 2 
U 


= 0 , 


where 


and 


v t = v f v 1 , /, 


vl 


3 | 3 i X — •> 

X 6 + C 6 vl u 


du 


X 


& 7 o ofv2i fv 2 1 o • 

dy K z y z 1 + xfv l 


ft 2 = c t3 exp (-c t4 i> 2 ) . 


By construction, these equations have the simple solution, 

u = Ku r y, 


S= — 

ny 


(26) 


(27) 

(28) 


(29) 


Transforming to wall units, y + = yu T /u , u + = u/u r: this solution becomes, 


X = M/+, S+ = 4-S=A ( 3 °) 

This solution is the extension of well known log-law behavior to the entire inner layer from the wall through 
the viscous sublayer and into the log-law region. In the log layer, Reynolds’ stresses dominate molecular 
stresses. With the introduction of the Boussinesq approximation and the log-law velocity distribution, both 
velocity gradient and eddy viscosity can be determined, 



du + 1 
dy+ Ky + ’ 


y + » l. 


(31) 


In developing the near-wall or viscous sublayer components of SA, this simple behavior was retained for the 
new solution variable u and the modified vorticity S by introducing the eddy viscosity correlation function 
/ v i, the definition of modified vorticity (via f v 2 ) and the re-definition of r (not shown). Note that the presence 
of the laminar suppression term, f t 2 , in SA is passive with respect to the simple solution; contributions from 
production and wall destruction cancel in (26). 

Substituting the simple solution (30) into either ^-momentum or the definition for S then gives, 


du + _ (? vX + (/q/ + ) 3 
dy + + (m/+) 3 ( 1 + k y+) 


( 32 ) 
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This equation can be integrated (via Mathematica 53 ) and the constant of integration determined from the 
no-slip boundary condition i/(0) = 0. Using complex arithmetic, the solution is, 


« + (A) = X) 

2=1 


Kzf( 3 + 4 Zi) 


[log ( ny + - Zi) - log (-Zi)] , 


where z, are the four solutions to the quartic equation, 

4i + z f + z t = o 


(33) 


(34) 


The solution can be simplified and rewritten using real arithmetic, 


u + (y + ) = B + ci log (( y + + a^ 2 + b\) - c 2 log (( y + + a 2 ) 2 + bf) 

- c 3 ArcTan[y + + ax, b\] - c 4 ArcTan [y + + a 2 , b 2 \, (35) 

where ArcTan[x, y\ is the Mathematica function equivalent to the Fortran function atan2(i/, x). For the 
values of ft = 0.41 and c v i =7.1, the constants are given by, 


B = 5.0333908790505579 
ai = 8.148221580024245 
a 2 = -6.9287093849022945 
ci = 2.5496773539754747 
c 3 = 3.599459109332379 


h = 7.4600876082527945 
b 2 = 7.468145790401841 
c 2 = 1.3301651588535228 
c 4 = 3.6397531868684494 


Analysis of this solution for large y + reveals that SA asymptotically produces a log-law with a shifted origin 
compared to the conventional formulas, 

du 1 

U + log ( y + + 1 /k) + B, — — as y + oo. (36) 

ft v y ay + ft^/ + + 1 

The shift in the asymptotic gradient can also be derived from (32). The origin shift is minor, but is easily 
noticeable in law of the wall velocity plots. 

Figure 24 shows the law of the wall velocity profile for SA (35) compared to its asymptotic form (36) and 
to Spalding’s composite formula, 


y + = u + + exp(— kB) 


Kp (ft^ + ) — 1 — KU + — -(ft^ + ) 2 — -(/Vft + ) 3 


(37) 


Here the values of ft = 0.41 and B = 5 have been plotted; this value of B is consistent with the calibration 
of SA (and in particular c v \ = 7.1). SA differs from Spalding by about 3.4% at y + = 50. 


u+ 
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Introduction 



Cartesian methods very useful for simulations with complex geometry 
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Introduction 


Viscous simulations remain a challenge due to meshing 
requirements near wall 

Primary approaches 
e Conformal layers near wall 
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Introduction 


Viscous simulations remain a challenge due to meshing 
requirements near wall 

Primary approaches 
0 Conformal layers near wall 

• Replace near-wall cells with body-fitted mesh 

• Sacrifice simplicity of pure Cartesian approach 



© Pure Cartesian Approaches 


Introduction 
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Introduction 


Viscous simulations remain a challenge due to meshing 
requirements near wall 

Primary approaches 

O Conformal layers near wall 



Pure Cartesian Approaches 

• Less widely explored 

• Contend with inefficiency of isotropically refined cells near wall 

• Confront mesh irregularity for 2nd-derivatives in Navier-Stokes 


RANS Solutions on Cartesian Cut-Cell Meshes 



Open issues/ToDo List 

• Accurate differencing for Navier-Stokes 

• Confront mesh irregularity at wall - smooth skin friction 

• Strengthen solver for 2nd-order PDEs 

• Implement initial turbulence model for RANS 

• Develop initial wall model for use within cut-cells 

• Demonstrate for simple cases 

• Examine limits (numerics or physics?) 
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Differencing for Viscous Flux 


Cut-Cell Discretization 


• Integration of Navier-Stokes requires 
accurate viscous fluxes on face-centroids 
of cut-cells 


• Viscous fluxes are assembled from the 
shear stresses 

r xx = 2\iu x — 2/3/xV • v 

7~xy Vx) 7~yx 

Tyy = 2[lVy ~ 2/3^V * V 



Differencing for Viscous Flux 



Cut-Cell Discretization 


• Mitigate mesh irregularity by recentering 
within cut-cells 


• Recenter from cell centroids (A,C) to per- 
pendicular through the face centroids 
(B,D) using gradients 




Use simple differencing on recentered 
data 


Ux 


Ub — Ud 

Ax 


• Accuracy study shows 2nd-order 
convergence for 2nd derivatives 
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Mesh Irregularity at Cut-Cells 



Confronting mesh irregularity 

• Literature shows well-known examples of poor/noisy skin friction without any clear 
path forward (Coirier’94, Marshall’02, Zhao’10) 


• Cf unacceptably irregular & shows poor mesh convergence 
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Mesh Irregularity at Cut-Cells 



Examine wall-normal gradients 

• Consider wall-normal component of gradient in non-aligned Cartesian cell 
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Mesh Irregularity at Cut-Cells 



Examine wall-normal gradients 

• Consider wall-normal component of gradient in non-aligned Cartesian cell 


• Imagine a function varying 
quadratically with distance from 
the wall... 

• Now, reconstruct this function with 
a linearly-exact gradient operator 
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Mesh Irregularity at Cut-Cells 



Examine wall-normal gradients 

• Consider wall-normal component of gradient in non-aligned Cartesian cell 

• Imagine a function varying 
quadratically with distance from 
the wall... 

• Now, reconstruct this function with 
a linearly-exact gradient operator 

• Reconstruction is noisy due to 
staggering of centroids 

• Gradient along wall will appear 
noisy, despite linearly exact 
reconstruction 



14 


Mesh Irregularity at Cut-Cells 


Regularizing the wall-flux 

The Basic Idea: 

• Recognize that cell-centroid data is correct, 
and perform quadratic reconstruction of the 
function through these values 

• Use analytic slope of the quadratic to evaluate 
wall-flux — do not use the cell-centroid slopes 



• Maintains full conservation 


• Exactly satisfies wall BC 



Viscous Cartesian Cut-Cells 

Regularizing the wall-flux 

• Quadratic reconstruction in wall-normal direction of cut-cells 

• Quadratics “regularize” the viscous flux stencil in the wall-normal direction 

• Turns out that need to “regularize” gradient stencils for 1-away cells too, (these 
cells have corrupted gradient stencils also) 


• Use quadratics anytime integrating-to-wall 



o2 



dU 


U 


RANS Equations 



Turbulence Model Implementation 

• Spalart-Allmaras turbulence model for jit 

• Fully turbulent (No trip terms) 

• “Negative SA Model” (S. Allmaras / T. Oliver, ’ 08 ) 

• Second-order advective terms in model 

• Details in paper 



• Sethian’s Fast Marching Method 
(FMM) for distance metric 

• Eikonal Eq. solver on 
multilevel cut-cell mesh 
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Cut-Cell Wall Model 



Use subgrid-scale wall modeling to mitigate cell-counts 





Body-fit with subgrid 
wall model 

• Craft (’04), 

• Myers & Walters (’05) 

• Bond & Blottner (‘07) 


Immersed Boundary 
with wall model 

• Capizzano (’10) 
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Cut-Cell Wall Model 



Use subgrid-scale wall modeling to mitigate cell-counts 


• Solve limiting form of RANS on subgrid 
in each wall-bounding cell 

• Doesn’t assume a form of solution, 
therefore more general than traditional 
(analytic) wall functions 

• Most approaches based on diffusion- 
model 

• Ultimate goal is develop PDE-based 
subgrid model, but for now, use 
analytic wall functions 
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Cut-Cell Wall Model 



Spalding Eq. gives implicit relationship for velocity u + , requiring iteration 

y + = u + + e~ kB (e ku+ - 1 - ku + - i( ku + ) 2 - i(£:tt + ) 3 ) 

• This relationship implies particular profile for eddy viscosity 

• May not exactly match exactly outer turbulence model (Kalitzin’05) 

• Mismatch can lead to issues with mesh convergence 

SA Wall Function: 

• Allmaras derived analytic relationship for u + (y + ) using limiting form of SA model 
for near-wall flow 

u + (y + ) = B + Cl log((j/ + + ai ) 2 + b\) - c 2 \og((y + + a 2 ) 2 + &|) 

c 3 arctan(^^-) - c 4 arctan(^_) 

• Explicit evaluation of u + , derivation in paper 

• Profiles consistent with outer flow (energy, momentum, eddy viscosity) 
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Cut-Cell Wall Model 



Use Spalart-Allmaras wall-function 


• Based on diffusion model of 
near-wall flow 

• Matches Spalding in sublayer 
and log-layer 

• Differs by 3.4% @ y + = 50 

• Matches wall-resolved SA 
calculations from wall to wake 



21 


Wall Model Implementation 


Wall model plays same 
role as quadratic did at 
low-cell Reynolds 
numbers 




• Each cut-cell gets a wall model describing variation in wall-normal 
direction 

• Wall model solved to match wall boundary-condition with solution 
on Cartesian mesh 

• Currently, use SA wall functions for wall model 
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Coupling: Cartesian -► Wall model 



Replace quadratic with 
wall model 



• Wall-functions have sensitivity to location of first point 


• Regularize the distance by choosing forcing points, F, a fixed 
distance h away, (following Capizzano) 

• Reconstruct for data at F from Cartesian mesh 


• Solve the ID wall model, (Newton iteration for friction velocity, u T ) 
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Coupling: Wall model -► Cartesian 



Replace quadratic with 
wall model 


Wall model provides 

• Velocities at centroid of cut-cell 

• Wall-normal gradients in cut-cells 

• Viscous wall flux in cut-cells 

• Tangential velocity components that feec- 
viscous flux calcs, at all the cut-faces 



D 


F 

© 

© 
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Results 



• Laminar flat plate 

• NACA 001 2 

• Turbulent Flat Plate 

• RAE 2822 

• Subsonic 

• Transonic x 2 
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Laminar Flat Plate Boundary Layer 



Mx = 0.5, Ret = 5000, a= 15 



Use quadratic reconstruction in wall-normal direction 
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Laminar Flat Plate Boundary Layer 



• Quadratic reconstruction in wall-normal direction 


• “Standard conclusions apply” 

• Need 12-20pts in boundary layer 

• Least dissipative flux function does best (e.g HLLC beats van Leer) 

• Limiters can clip the “viscous overshoot” but contribute to skin-friction 
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Laminar Flat Plate Boundary Layer 


Mao = 0.5, Ret = 5000, 
a= 15° 


No Limiter 


- 


- 


» — • Rc x = 5000 (-13 cells) 



■ 

■ — ■ Re x = 10000 (-17 cells) 


■ 

- 

4 — ♦ Re x = 50000 (-40 cells) 


- 





Laminar Flat Plate Boundary Layer 



Mx = 0.5, Rei - 5000, a - 15° 
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Laminar NACA 001 2 



Mao = 0.5, Rec = 5000, a= O' 



Quadratic reconstruction in the wall-normal direction, Ax « 0.1 6%c 
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Laminar NACA 001 2 



M» = 0.5, Rec - 5000, a= 0° 



• Quadratic reconstruction in the wall-normal direction, Ax « 0.1 6%c 

• Separation at x = 81 .6%c 
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Laminar NACA 001 2 



Mx = 0.5, Rec - 5000, a= 0° 


Multigrid performance 
with 5 coarse grids 
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Laminar NACA 001 2 



M» = 0.5, Rec - 5000, a= 0° 



• Smooth skin friction, separation at x = 81 .6%c 

• Good comparison with published data 
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Turbulent Flat Plate Boundary Layer 



Moo = 0.5, Rei = 5 x 1 0 5 , a = 0° 

• Integrate to the wall using quadratic wall-normal reconstruction (low 
cell-Reynolds number) 
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Turbulent Flat Plate Boundary Layer 



Moo = 0.5, Ret = 5 x 1 0 5 , a = 0° 


• Integrate to the wall using quadratic wall-normal reconstruction (low 
cell-Reynolds number) 

• Mesh converged solutions @13 refinements, Ax = 3.6 x 10' 4 , y + « 8.8 



Tangential Velocity 
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Turbulent Flat Plate Boundary Layer 



Moo = 0.5, Ret = 5 x 1 0 5 , a = 0° 


• Integrate to the wall using quadratic wall-normal reconstruction (low 
cell-Reynolds number) 

• Mesh converged solutions @13 refinements, Ax = 3.6 x 10' 4 , y + « 8.8 



36 


Turbulent Flat Plate Boundary Layer 



Moo = 0.5, Ret = 5 x 1 0 5 , a = 0° 


Examine performance of wall model 
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Turbulent Flat Plate Boundary Layer 



Moo = 0.5, Ret = 5 x 1 0 5 , a = 0° 


Examine performance of wall model 
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Turbulent Flat Plate Boundary Layer 



Moo = 0.5, Re L = 5 x 1 0 5 , a- 15° 


Compare with wall-aligned results 
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Turbulent Flat Plate Boundary Layer 



Moo = 0.5, Re L = 5 x 1 0 5 , a- 15° 


Compare with wall-aligned results 
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RAE 2822 - Resolution Requirements 


• Subsonic & transonic examples 


• All used SA wall-model 

• Meshes 

• Coarse: Ax = 0.1 %c 

► y + « 200-350, ~50k cells 

• Medium: Ax = 0.05%c 

► y + ~ ] 00-1 70, —11 5k cells 

• Fine: Ax = 0.025%c 

► y + « 50-80, ~150k cells 




RAE 2822 - AGARD Case 1 



Moo = 0.676, Re c = 5.7 x 1 0 6 , 


• Experimental a = 2.4° 

• Simulation a = 0.1 .85° 

• Used SA wall-model 

• No circulation correction 


C L = 0.566° 
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RAE 2822 - AGARD Case 1 


Moo = 0.676, Re c = 5.7x1 0 6 , Cl = 0.566' 


Multigrid performance 

• Convergence of medium grid 

• 4 levels of grids 

• Typical of performance on 
RAE examples 
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RAE 2822 - AGARD Case 1 



Moo = 0.676, Re c = 5.7 x 1 0 6 , C L = 0.566° 


• Slight differences between 
grids near suction peak 

• Results largely invariant 
with mesh 

• Agreement comparable to 
other simulations in 
literature 
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RAE 2822 - AGARD Case 1 



M» = 0.676, Re c = 5.7 X 1 0 6 , Cl = 0.566° 


• Much more variation 

• Similar C/s near T.E where 
B-L is thick 

• Results degrade moving 
forward 

• Near L.E. B-L is thin, only 
fine mesh predicts data 

• Consistent with y + 
requirements of wall model 
in flat plate examples 
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RAE 2822 - AGARD Cases 6 & 1 0 


AGARD 
Case 6 




Mo = 0.725, Re c = 6.5 x 1 0 6 , C L = 0.743° 


M„ = 0.75, Re c = 6.2 x 1 0 6 , C L = 0.743° 


RAE 2822 - AGARD Cases 6 & 1 0 



AGARD 
Case 6 



Mo = 0.725, Re c = 6.5 x 1 0 6 , C L = 0.743° 


AGARD Case 6 

• Experimental a = 2.92° 

• Simulation a = 2.44° 

• Shock @ ~0.5c 

• No circulation correction 
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RAE 2822 - AGARD Cases 6 & 1 0 



AGARD Case 1 0 

• Experimental a = 3.19° 

• Simulation a = 2.54° 

• Stronger shock @ ~0.58c 

• No circulation correction 
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RAE 2822 - AGARD Cases 6 & 1 0 





Mco = 0.725, Re c = 6.5 x 1 0 6 , Cl = 0.743° 


Mo o = 0.75, Re c = 6.2 x 1 0 6 , Cl = 0.743° 
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RAE 2822 - AGARD Cases 6 & 1 0 





Mco = 0.725, Re c = 6.5 x 1 0 6 , Cl = 0.743° 


Mo = 0.75, Re c = 6.2 x 1 0 6 , Cl = 0.743° 
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RAE 2822 - AGARD Cases 6 & 1 0 




Mco = 0.725, Re c = 6.5x1 0 6 , Cl = 0.743° 


Mo o = 0.75, = 6.2 x 1 0 6 , Cl = 0.743° 
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RAE 2822 - AGARD Cases 6 & 1 0 





• Cp distributions largely invariant with mesh 


• Agreement comparable with other published solutions 

• Wall spacing on coarse mesh is 0.1 %c 
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RAE 2822 - AGARD Cases 6 & 1 0 




• Close-up of shock-boundary layer interaction 
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RAE 2822 - AGARD Cases 6 & 1 C 


AGARD 
Case 6 



Mach Contours 





Close-up of shock-boundary layer interaction 
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RAE 2822 - AGARD Cases 6 & 1 0 




Moo = 0.725, Re c = 6.5 x 1 0 6 , C L = 0.743° M„ = 0.75, Re c = 6.2 x 10 6 , C L = 0.743° 
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Summary 



Progress 

/ Develop accurate differencing for Navier-Stokes 

/ Confront mesh irregularity at wall for smooth skin friction 

/ Robust multigrid convergence for RANS Eqs. 

/ Implemented Negative SA turbulence model 

/ Develop initial wall model for use within cut-cells - SA wall function 

/ Demonstrate for simple cases, including shock B-L interaction 

/ Resolution requirements currently driven by analytic wall model, 
achieve Cp with Ax « 0.1 %c, but C/ needs y+ in log-layer 

® Immediate goal is to develop PDE-based wall model with more 
extensive physics 
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